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Abstract 

The commensal microbiota offish skin is suspected to provide a protection against opportunist infections. The skin offish 
harbors a complex and diverse microbiota that closely interacts with the surrounding water microbial communities. Up to 
now there is no clear evidence as to whether the host regulates the recruitment of environmental bacteria to build a specific 
skin microbiota. To address this question, we detected Quantitative Trait Loci (QTL) associated with the abundance of 
specific skin microbiota bacterial strains in brook charr {Salvelinus fontinalis), combining 16S RNA tagged-amplicon 454 
pyrosequencing with genetic linkage analysis. Skin microbiota analysis revealed high inter-individual variation among 86 F2 
fish progeny based upon the relative abundance of bacterial operational taxonomic units (OTUs). Out of those OTUs, the 
pathogenic strain Flavobacterium psychrophilum and the non-pathogenic strain Methylobacterium rhodesianum explained 
the majority of inter-individual distances. Furthermore, a strong negative correlation was found between Flavobacterium 
and Methylobacterium, suggesting a mutually competitive relationship. Finally, after considering a total of 266 markers, 
genetic linkage analysis highlighted three major QTL associated with the abundance of Lysobacter, Rheinheimera and 
Methylobacterium. All these three genera are known for their beneficial antibacterial activity. Overall, our results provide 
evidence that host genotype may regulate the abundance of specific genera among their surface microbiota. They also 
indicate that Lysobacter, Rheinheimera and Methylobacterium are potentially important genera in providing protection 
against pathogens. 
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Introduction 

The study of the beneficial effects of bacteria and their influence 
on host health is a growing field. Namely, research that has 
explored host microbiota variability in space and time suggests the 
presence of a host genetic component into the development of 
mutualism communities [1,2]. Similar conclusions were found in 
several studies on twins and on their core gut microbiota [3,4,5,6]. 
Furthermore, microbiome homeostasis seems to be the key to 
resistance against some diseases previously considered exclusively 
influenced by genetic factors [7]. To determine precisely which 
genes are involved in the recruitment of specific bacterial strains, 
some studies looked at gene expression in presence of symbiotic 
bacteria. Recent studies specifically targeted genes already known 
as being involved in innate or adaptive immunity (for example 
IgA) [8,9,10]. 

In fish, the influence of host genetic background on commensal 
bacterial community structure is poorly known. However, 
advances in the field of probiotic development indicate that 
endogenous bacteria are able to outcompete pathogens [11,12]. 
Results obtained in aquaculture settings are consistent with 
previous work showing that non-pathogenic bacteria associated 



with animal mucosa can contribute to the host health by providing 
protection against pathogen infections [13,14,15,16]. The genetic 
basis of the host immune response, especially at the major 
histocompatibility complex (MHC), toll-like receptor (TLR) and 
immunoglobulin loci, has been well studied [17,18,19,20]. 
However, aside from the immune response, the influence of the 
host on the development of bacterial symbiosis remains poorly 
understood. 

Here, the main objective of the present work was to investigate 
the potential link between host genotype and skin microbiota 
composition. As a host model, we focused on brook charr 
(Salvelinus fontinalis), a salmonid that harbors a complex and 
dynamic skin microbiota [21,22]. The first specific objective 
sought to document the genetic variation present in the host 
population that might underpin the inter-individual variations in 
the skin microbiota. The second specific objective was to 
characterize the relationship (cooperation/competition) existing 
among the different bacterial Operational taxonomic units (OTU) 
within the skin mucus, which influences the overall structure of the 
skin microbial community. The third objective was to identify host 
genetic regions associated to the abundance of specific skin 
bacterial OTUs (i. e. quantitative trait loci (QTL) associated with 
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Table 1. Descriptive statistics of sequencing to determine the microbiotoa for each individual. 





Sample ID 


number of reads 


coverage 


npshannon index 


Number of 
OTU total 


Number of OTU >10 
reads 


92 


229 


0.955975 


1 .690045 


76 


19 


93 


1584 


0.99131 


0.436879 


74 


28 


94 


669 


0.977401 


1.507996 


242 


20 


95 


738 


0.986166 


1.416086 


226 


27 


96 


503 


0.953079 


1.395581 


149 


29 


97 


385 


0.946203 


1.774329 


86 


33 


98 


955 


0.98389 


0.663968 


79 


30 


99 


770 


0.978852 


1 .306709 


114 


29 


100 


1018 


0.985946 


0.658787 


107 


30 


101 


1482 


0.994792 


0.842304 


188 


21 


102 


1667 


0.992283 


0.691189 


96 


31 


103 


711 


0.99179 


0.946539 


97 


23 


104 


1072 


0.987599 


0.719994 


159 


25 


105 


1662 


0.991645 


0.371924 


73 


26 


106 


1304 


0.995734 


0.658557 


93 


18 


107 


1617 


0.99416 


0.290016 


49 


21 


108 


1237 


0.990451 


0.634742 


73 


30 


109 


1063 


0.987778 


0.679124 


137 


26 


110 


752 


0.97281 


0.656197 


105 


29 


111 


731 


0.986217 


0.63276 


79 


22 


112 


628 


0.975779 


0.630137 


63 


23 


113 


892 


0.983549 


0.446389 


49 


24 


114 


1333 


0.993127 


0.713271 


127 


22 


115 


1109 


0.983313 


1 .020009 


266 


35 


116 


532 


0.96761 1 


2.709733 


269 


32 


117 


649 


0.963145 


2.108945 


237 


41 


118 


404 


0.923977 


2.645631 


222 


29 


119 


969 


0.97625 


1.019139 


177 


40 


120 


581 


0.967647 


1.975631 


205 


26 


121 


468 


0.948413 


2.692236 


202 


34 


122 


744 


0.980551 


0.876389 


117 


27 


123 


953 


0.983204 


0.832745 


133 


32 


124 


530 


0.967846 


2.29117 


186 


30 


125 


797 


0.977901 


0.703269 


79 


30 


126 


411 


0.97561 


0.904985 


116 


18 


127 


242 


0.76 


3.452235 


182 


32 


128 


1230 


0.989637 


0.504536 


69 


26 


129 


831 


0.978723 


1.537117 


115 


41 


130 


1449 


0.983421 


0.548428 


123 


38 


131 


499 


0.969697 


2.118408 


144 


29 


132 


431 


0.945813 


2.133064 


181 


25 


133 


494 


0.96114 


1 .408929 


117 


35 


134 


746 


0.980254 


0.508333 


38 


22 


135 


662 


0.977887 


2.41624 


180 


36 


136 


754 


0.987755 


0.956343 


118 


18 


137 


2008 


0.993154 


0.557917 


74 


30 


138 


714 


0.976619 


1.632525 


101 


36 


139 


798 


0.967632 


1.534968 


129 


35 


140 


469 


0.917197 


3.267019 


261 


35 
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Table 1. Cont. 







Sample ID 


number of reads 


coverage 


npshannon index 


Number of 
OTU total 


Number of OTU >10 
reads 


141 


224 


0 980897 


0 667835 


33 


5 


142 


517 


0 957377 


7 388776 


1 80 


36 


143 


483 


0 9471 95 


2 44671 1 


1 57 


36 


1 44 


414 


n Q^RRAR 


3 1 24298 


1 88 


38 


1 45 


930 


n Q7^AA3 


1 391 363 


1 68 


43 


146 


769 


n Q7Qfifi1 


7 07A70R 


1 40 


28 


1 47 


1219 


0 9951 22 


1 fiAIARQ 

1 .OH- 1 to? 


118 


-| 7 


148 


1 254 


n QR7AA 


1 127563 


1 48 


40 


1 49 


1 074 


n Q77ft7R 


7 ^nonziQ 


318 


55 


1 50 


1 108 




2 793443 


321 


47 


1 51 


3082 


0 9Q8643 


0 332436 


58 


1 6 


1 52 


1429 


0 9Q6561 


1 R70656 


203 


25 


1 53 


1 955 


0 QCK7Rfi 


0 81 5557 


1 67 


27 


1 54 


1 449 


U. jOjjDD 


1 74373 


354 


54 


1 55 


1 074 


0 97361 1 


1 716763 


321 


41 


1 56 


1016 


0 QRQ31 ft 

U.703J 1 D 


0 ^RR71 Q 

U.JOO/ I 3 


62 


21 


1 57 


1 483 


0 QRAR^fi 


1 .J / HUUJ 


216 


46 


1 58 


2004 


fl QQRQAA 


0 791 1 57 


44 


1 0 


1 59 


2043 


0 997359 


0 R1 7388 

U.O 1 ZJOO 


92 


1 8 


160 


2468 


0 995741 


0 305525 


70 


23 


161 


1 215 


0 97Q566 


1 928985 


346 


48 


162 


1 904 


U.77JZOO 


n R71 ftOQ 

U.OZ 1 DU3 


1 64 


37 


1 63 


1 046 


0 985799 


1 569251 


1 88 


40 


164 


1 1 45 


0 QQ3R1 A 


n 7AA c ;r7 

U./ HtJO/ 


98 


1 9 


1 65 


995 


0 QR1 AR1 


1 A7 1 1 0A 

I .H-Z. I I UH- 


243 


37 


166 


725 


0 97921 


2 459141 


211 


36 


1 67 


727 


0 QRR7A^ 


2 142562 


1 3 1 


24 


168 


841 


0 9RQ0R 


2 297402 


1 51 


31 


169 


552 


0 9681 16 


2 76388 


1 78 


35 


170 


1204 


0.99177 


1.912832 


162 


32 


172 


2010 


0.993946 


0.610072 


155 


36 


173 


1454 


0.987296 


1.893132 


316 


51 


174 


1056 


0.979651 


2.949199 


335 


53 


175 


824 


0.95098 


2.845177 


368 


46 


176 


851 


0.96729 


2.039242 


202 


47 


177 


778 


0.985507 


2.23077 


221 


25 


178 


2142 


0.998093 


0.359619 


52 


13 



The coverage was estimated by a Good's estimator index, npshannon: non parametric index of Shannon. OTU classification was done with a treshold of 97% identity. 
doi:1 0.1 371 /journal.pone.01 02649.t001 



the abundance of each constituent bacterial genus). Accordingly, 
we analyzed the taxonomic structure of skin bacterial community 
using 16S tagged-amplicon pyrosequencing in 86 F2 fish progenies 
obtained in the previous work of Sauvage el al. (2012). Inter- 
individual variation in abundance of each bacterial strain detected 
in host microbiota was afterwards projected on the results of the 
genetic linkage map in order to identify quantitative trait loci 
(QTL) of specific skin bacterial OTUs. 



Materials and Methods 

Ethics statement 

All fish were reared and the experiment conducted stricdy 
following guidelines required by the "Comite de Protection des 
Animaux de FUniversite Laval (CPAUL, http://www.vrr. ulaval. 
ca/deontologie/cpa/index.html?accueil). The CPAUL reviewed 
and approved all experimental procedures used in this study. 
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Class 



Otu 




Figure 1. Taxonomic structure of the bacterial community of 
fish skin mucus at three different taxonomic levels: Phylum, 
Classl and OTU. a) Proteobacteria; b) Bacteroidetes, c) Alphaproteo- 
bacteria, d) Flavobacteria, e) Gammaproteobacteria, f) OTU 50 
{Methylobacterium rhodesianum), g) OTU 36 (Flavobacterium psychro- 
philum). 

doi:1 0.1 371 /journal.pone.01 02649.g001 



Fish sampling 

The study population was generated using two different 
populations of brook charr. The first one (D) is a domestic 
population used in aquaculture for more than 100 years 



(Pisciculture de la Jacques Cartier - Cap-Sante, Quebec, Canada). 
The other population (L), is an anadromous population from the 
Laval River near Forestville (North of the St. Lawrence River, 
Quebec, Canada). Breeders from the L population were kept in 
captivity at the Station aquicole de 1'ISMER (Quebec, Canada, 
48°31'N, 68°28'W) under natural photoperiod and temperature. 
Crosses between 10 sires of each population (L and D; F 0 
generation) with 10 dams (L and D) were performed to generate 
10 full-sib outbred hybrid families (LxD - F) generation). 
Subsequently, six Fj individuals were bi-parentally crossed to 
obtain three F 2 families. The F 2 family selected for the present 
study demonstrated the lowest post-hatch mortality rate (<2%) 
and the greatest abundance. Fish were raised in the same tank and 
had fasted for 12 hours before sampling. Fish were measured 
(mean = 28.8± 1.77 cm) and weighed (mean = 276.7±59.48 g) to 
calculate Fulton index [23]. 

Mucus was sampled using a sterile swab on the same area for all 
the fish [24] . We choose to sample the skin between the adipose fin 
and the caudal fin because this area was undisturbed by fish 
handling. Samples were stored in a sterile micro-centrifuge tube 
containing lysis buffer (Tris 50 mM, EDTA 40 mM, sucrose 
0.75 g) and immediately stored in — 80°C until DNA extraction. 

DNA extraction 

DNA was extracted using a modified protocol of salt-extraction 
from Aljanabi & Martinez (1997). During the first lysis step, 
22.6 |iL of lysozyme (40 mg/mL) was added to the sample and 
incubated 45 minutes at 37°C. After this step, 22.6 [lL of 
proteinase K (20 mg/uL) and 90 uL of SDS (10%) were added 
to the lysate and incubated at 55°C over night with agitation. The 
aqueous phase was transferred into a clean Eppendorf tube 
containing 600 uL of NaCl 6M, mixed and centrifuged 20 min at 
14,800 rpm. The supernatant was transferred again into a clean 
Eppendorf tube containing 1 volume of cold isopropanol, mixed 
and stored 30 minutes at — 20°C. The mixture was centrifuged 20 
minutes at 14,800 rpm and the supernatant was thrown away. 
The pellet was washed with cold ethanol 70%, air-dried and finally 
resuspended in 25 |iL of sterile MilliQ H 2 0. Subsequendy, DNA 
integrity and quantity were controlled using a Nanodrop 
instrument (ND-1000, Nanodrop). 

Microbial 16S pyrosequencing 

Each DNA sample, the 16S gene was PCR amplified using 
Takara Ex taq premix (Fisher). All PCR reactions were performed 
in a reaction volume of 50 |jL containing 25 |iL of premix Taq, 
1 |xM of each primer and sterile MilliQ H z O to up to 50 uL. A 
general reverse primer (R519) combined with B primer (Roche) 
was used for amplifications in combination with one of 86 
uniquely tagged forward primers (F63-targeted) combined with A 
primer (Roche) (for primer sequences see [25,26]. The mean 
length of the amplified fragment was 450 bp. This procedure 
facilitates the parallel sequencing of thousands of samples on the 
same run and to reassign each reads to their respective samples. 
PCR conditions were applied as follows: after a denaturing step of 
30 s at 98°C, samples were processed through 30 cycles consisting 
of 10 sec at 98°C, 30 sec at 55°C, and 30 sec at 72°C. The final 
extension step was done at 72°C for 4 min 30 sec. Following the 
amplification step, samples were purified using AMPure Beads 
(Beckman Coulter Genomics). Samples were adjusted to 100 uL 
with EB (Qiagen), 63 uL of beads were added. Samples were 
mixed and incubated for 5 min at RT. Using a Magnetic Particle 
Concentrator (MPC), the beads were pelleted against the wall of 
the tube and supernatant was removed. The beads were washed 
twice with 500 uL of 70% ethanol and incubated for 30 sec each 
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Table 2. Taxonomic variation observed between male and female. 





Condition 


N° OTU 


Best Hit on Blast 


Identity 


Relative abundance 


More abundant in female 


82 


Amaricoccus kaplicensis 


98% 


0.00037647 




115 


Loktanella salsilacus 


99% 


0.00015338 




6 


Methylovirgula ligni 


91% 


0.00026492 




163 


Micrococcus antarcticus 


99% 


0.00061351 




88 


Paracoccus yeei 


94% 


0.00019521 




75 


Pseudochrobactrum kiredjianiae 


99% 


0.00050196 




128 


Psychrobacter faecalis 


99% 


0.00100392 




131 


Rheinheimera pacifica 


99% 


0.00057168 




7, 11 


Singularimonas variicoloris 


89% 


0.00037647 




104 


Sphingopyxis witflariensis 


99% 


0.00015338 




113 


Stella humosa 


1 00% 


0.00018126 




142 


Stenotrophomonas nitritireducens 


99% 


0.00022309 


More abundant in male 


186 


Flavobacterium aquatile 


97% 


0.00029281 




190 


Hyphomicrobium sulfonivorans 


95% 


0.00039041 




4 


Methylosinus sporium 


92% 


0.00027887 




176 


Polynucleobacter necessarius subsp. Asymbioticus 


99% 


0.00015338 




129 


Pseudomonas hibiscicola 


99% 


0.00075294 




12 


Pseudorhodoferax soli 


85% 


0.00018126 




140 


Psychrobacter alimentarius 


99% 


0.00015338 




116 


Rhodobacter blasticus 


95% 


0.00016732 




175 


Simplicispira limi 


99% 


0.00034858 



This table summarize the OTUs differentially abundant between skin mucus communities between male and female (Metastat using a FDR correction, p-value<0.01). 
doi:1 0.1 371 /journal.pone.01 02649.W02 



time. Supernatant was removed and beads were air dried for 
5 min. Tubes were removed from the MPC and 24 |jL of EB were 
added. Samples were vortexed in order to suspend the beads. 
Finally, using the MPC, the beads were pelleted against the wall 
once more and supernatant was transferred to a new clean tube. 
Samples were quantified with Nanodrop before the amplification 
step. Amplicons were then quantified with Quant-iT PicoGreen 
dsDNA Assay Kit and mixed equally before being sent to the 
Plateforme d'Analyses Biomoleculaires (Institut de Biologie Inte- 
grative et des System.es, Universite Laval) for pyrosequencing on a 
454 GS-FLX DNA Sequencer with the Titanium Chemistry 
(Roche), according to manufacturer's procedure. 

16S sequence analysis 

All sequences are available on MG-RAST server (MG-RAST 
IDs: 4539915.3). The data were analyzed in two steps. First, CLC 
Genomics Workbench 3.1 (CLC Bio, Aarhus, Denmark CLC 
workbench BIO) was used to trim sequences for quality, recover 
and remove the primers' sequences and tags (minimum average 
quality score: 35 for a window of 50, number of differences to the 
primer sequence = 0, maximum number of differences to the 
barcode sequence = 0, number of ambiguous base calls = 0, 
maximum homopolymer length =8). In a second step, pre- 
processing and analysis were performed using the MOTHUR 
software [27]. All datasets were checked for chimeras with the 
chimera slayer algorithm implemented in MOTHUR. Standard- 
ization of the different samples was done by using the zscore which 
calculates the normalized abundance as follows: normalized 
abundance = (abundance - mean) / standard deviation 



[3,28,29]. An alternative method for normalization was also 
tested, which consisted of subsampling equal number of reads from 
each sample [30]. This method greatly reduced the numbers of 
sequences (224 sequences per samples), but gave the same results 
as zscore (Figure S2, Figure S3, Table SI). Therefore we preferred 
to keep zscore normalization to base our conclusions on a larger 
dataset. We used the Operational Taxonomic Unit-based method 
described by Costello el al. (2009) because it is not biased towards 
a predefined taxonomy. One index was retained to assess the 
quality of pyrosequencing: the sequence coverage index (Good's 
estimator). The sequence coverage index is a metric used to 
estimate the quality of the depth sequencing. All sequences were 
clustered into Operational Taxonomic Unit (OTU) using a 97% 
identity threshold and OTU were classified from phylum to genus 
using the program MOTHUR with the default settings. For all 
interesting OTUs (explaining the inter-individual variation or 
linked to a genetic region), we extracted all the sequences classified 
in that OTU and used the single best BLAST hit to identify the 
different species. Statistical differences in the abundance of each 
OTU were calculated with the software for metagenomic analysis 
(Metastats). 

To visualize potential differences across host progenies in terms 
of mucus bacterial community structure, distances between those 
communities were computed using the Yue & Clayton measure of 
dissimilarity (Thetayc). This index developed by Yue and Clayton 
(2005) is a measure of dissimilarity between the structures of two 
communities, meaning that this calculator takes in account the 
abundances of each OTU. Then, Distances were represented 
using a dendrograms based on this index and statistical robustness 
of the dendrogram was determined by a Unifrac weighted test. 
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Group 1 



Group 2 



Group 3 



Figure 2. Dendrogram analysis based upon ThetaYC index of 
bacteria found on the skin of the 86 brook charr individuals. 

Groups are defined with the Weighted Unifrac distance. The first and 
second groups are composed of closely related bacterial communities. 
The third group is an assemblage of dissimilar communities and is 
considered as an outgroup. The most differentiated groups are groups 
1 and 3 (Unifrac Score: 0.710811, p<0.0010) followed by the distance 
between groups 2 and 3 (Unifrac Score: 0.685361, p<0.0010), and the 
distance between groups 1 and 2 (Unifrac Score: 0.401674, p<0.0010). 
doi:1 0.1 371/journal.pone.01 02649.g002 

This test allows determining whether any of the samples have a 
significantly different structure than the other samples. A random 
(Monte Carlo) permutation test was performed to test whether or 
not the distance between two communities was greater than 
expected by chance alone. 

To visualize the distances between communities, we also applied 
a Principal Coordinate Analysis (PCoA) based on the Yue & 
Clayton measure of dissimilarity (Thetayc) using an eigenvector- 
based approach to represent multidimensional data in as few 
dimensions as possible. This method allows determining which 
OTUs are responsible for shifting the samples along the PCoA 
axes by measuring the correlation of the relative abundance of 
each OTU clustered in genera with the axes of the PCoA dataset. 
Statistical differences in the abundance of each OTU was 
calculated to highlight which genera were the causative agents of 
the differentiation between groups obtained in the dendrogram 
based upon the TefhaYC index [31]. All statistical analyses and 
graphics were carried out in the R environment (http:/ /www.r- 
project.org). 

QTL detection 

QTL analysis was carried out using the [R] package J/qtl [v. 
1.3.1, August 2012, http://churchill.jax.org/software/jqtl.shtrnl/]. 
QTL were projected on the consensus brook charr linkage map (see 
Sauvage et al., 2012) according to the following steps. (1) A single 
QTL analysis was performed using the Haley-Knott (HK) 
regression method (10,000 permutations) to reveal which LGs were 
carrying QTL. The most probable position of the QTL was defined 
at the position giving the largest logarithm of odds (LOD) score 
indicating the QTL was fixed. (2) A two QTL model based upon the 
Haley-Knott regression was used to refine the QTL detection across 
the genome with a resolution of 5 cM and eventually to detect two 
QTL on a single LG linked to a particular trait. (3) The best model 
fitting our data was used to compute the percent variance explained 
(PVE expressed in %) by the QTL. The chromosome-wide and the 
genome-wide thresholds were calculated for each LG using 10,000 
permutations. The 1.5 LOD confidence intervals were determined 
for all analyses following the Bayesian method implemented in the 
"bayesint" function in R/qtl. The bayesfnt function calculates an 
approximate interval (end points around the maximum LOD) for a 
given chromosome using the genome scan output. 

Results 

Mucus samples were obtained from 86 F2 fish progenies (44 
males and 42 females) issued from the same family. A total of 
87,940 high-quality, classifiable sequences were obtained with an 
average number of 1022 ±540 per sample, which were subse- 
quently classified in OTUs [1]. All of these sequences were 
successfully clustered into 9520 OTUs with 97% identity. A final 
trimming step was performed to focus on the most abundant OTU 
(which are represented by at least 10 sequences for the whole 
project) and result in a dataset containing 71,719 sequences (81% 
of the initial data set) clustered in 192 OTUs. The depth of 
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Bacterial diversity 




Group 

Figure 3. Alpha diversity of each group. Each group was defined by the dendrogram analysis based upon ThetaYC index (see figure 2). Alpha- 
diversity was calculated by the non-parametric index of Shannon. 
doi:1 0.1 371 /journal.pone.01 02649.g003 



sequencing and the coverage were estimated using the Good's 
estimator index; the coverage found was always higher than 90% 
except for one sample (#127), which exhibited 76% of coverage 
(Table 1, Figure SI). The number of OTUs per sample was not 
equally distributed and the mean number was 1 56.5 ±8 1 .7 
(ranging from 33 to 368). The alpha-diversity was estimated by 
calculating the non-parametric index of Shannon (npShannon) 
because of the non-equally distributed abundance of each OTU in 
the samples. The npShannon index ranged from 0.29 to 3.45, with 
a mean index of 1.43±0.83. 

Phylum, class, and OTU composition of all samples were 
represented in the Figure 1. Most of the OTUs belong to the 
Proteobacteria phylum (88.7%) and the Bacteroidetes phyla (9.7%). 



At the class level, most of the OTUs were classified as Alpha- 
proleobacteria (78.9%), Gamma-proteobacteria (6.5%) and Flavo- 
bacteria (13%). At the OTU level, the most abundant sequence 
was identified via BLAST as M. rhodesianum (69.8%) followed by 
F. psychrophilum (8%) (Figure 1 and Table SI). 

Male and female brook charr harbored slighdy different 
microbial communities. The Unifrac weighted distance between 
samples from males and females was statistically significant 
indicating that structures of bacterial communities were differen- 
tiated according to the sex of the host (Unifrac Score: 0.241695; p- 
value <0.001). However, there was not a sex difference in the two 
most abundant species, as this difference was based on 22 OTUs 
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Group 1 




Figure 4. Descriptive analysis of the Metastats. Numbers indicates the numbers of OTUs differentially present in the different groups of 
individuals. For details on the identities of the OTUs, please refer to the table 3; 4 and 5. 
doi:1 0.1 371 /journal.pone.01 02649.g004 



classified on 20 genera and 21 species, all of them being less 
abundant than 0.1% (Table 2). 

In terms of community structure, clustering based on ThetaYC 
index clearly defined two groups of samples, whilst a third group 
was less well defined and composed of highly variable communities 
(Figure 2). There was little variation among samples in the first 
and second group (group 1 and group 2). The third group (group 
3) was an assemblage of dissimilar communities and is considered 
as an external group. The variation of the alpha-diversity 
visualized by the npShannon index (Figure 3) indicated that 
alpha-diversity increases from group 1 to group 3. Furthermore, 
the bacterial communities within the third group were more 
diverse than groups 1 and 2. The same clustering was found with 
the normalization by subsampling (Figure S2, Figure S3, table SI). 
The difference between each group was greater than the difference 
between males and females (Figure 2). The most differentiated 
groups were groups 1 and 3 (Unifrac Score: 0.710811, p-value < 
0.001) followed by the distance between groups 2 and 3 (Unifrac 
Score: 0.685361, p-value <0.001), and the distance between 
groups 1 and 2 (Unifrac Score: 0.401674, p-value <0.001). This 



was correlated with the numbers of OTU that were differentially 
abundant between those 3 groups (Figure 4). The divergence 
between groups 1 and 2 was based upon 54 OTU classified in 26 
genera and 41 species (Table 3). The microbiota of the first group 
was dominated by M. rhodesianum, and host individuals from this 
group were tighdy clustered. The second group was more diverse, 
dominated by Flavobacterium and had a lower abundance of M. 
rhodesianum. For groups 1 and 3, the difference was based upon 
109 OTU classified in 45 genera and 70 species (Table 4). Finally, 
26 OTU were differently abundant between group 2 and 3 and 
classified in 17 genera and 22 species (Table 5). 

PCoA analysis displayed the same patterns of divergence 
between the three groups (Figure 5). Each group found in the 
Unifrac analysis was highlighted in the PCoA. To understand 
which OTUs were responsible for the differentiation of the cluster 
on the PCoA axes 1 and 2, individual correlation coefficients were 
calculated. Three OTUs were highly correlated with the two axes 
(a>0.6); OTU 36, 50 and 149. Those OTU belonged to two 
species: Flavobacterium psychrophilum (OTU 36) and M. 
Rhodesianum (OTU 50 and 149). The negative correlation 
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Table 3. Detailed information related to the OTU differentially abundant between individuals belonging to group 1 and group 2. 



N° OTU 


mean sbundBncc in group 1 


1116311 abundance in group 2 


dcsi nil on Did ST. 


36 


0.00484 


0.037926 


Flavobocterium psychrophilum 


122 


0.000778 


0.010066 


Pseudomonas cedrina subsp. Fulgido 


37 


0.000471 


0.006593 


Flavobocterium terrigena 


68 


0.0014 


0.006629 


Methylocelta palustris 


128 


0 


0.004734 


Psychrobacter faecaiis 


69 


0.000531 


0.004942 


Pseudorhodobacter ferrugineus 


190 


0 


0.004394 


Hyphomicrobium sulfonivorons 


139 


0.000153 


0.003833 


Lysobocter capsici 


70 


0.000831 


0.0043 


Methylobacterium salsuginis 


29 


0.000028 


0.003344 


Methylobacterium odhaesivum 


67 


0.000518 


0.003772 


Rhodobacter blasticus 


127 


0.000218 


0.003382 


Acinetobacter haemolyticus 


21 


0.000117 


0.002977 


Flavobocterium aquatile 


74 


0 


0.002743 


Paracoccus yeei 


52 


0.000331 


0.00291 1 


Sphingopyxis aloskensis 


126 


0.000969 


0.003524 


Pseudomonas xanthomarina 


89 


0 


0.002528 


Sphingomonas dokdonensis 


27 


0.000113 


0.001814 


Sphingopyxis chilensis 


63 


0.00071 


0.002374 


Paracoccus haeundaensis 


141 


0.000156 


0.001642 


Acinetobacter johnsonii 


84 


0.000079 


0.001484 


Pseudorhodobacter ferrugineus 


38 


0.000063 


0.001334 


Sejongia jeonii 


138 


0.000014 


0.001205 


Acinetobacter junii 


192 


0 


0.001099 


Methylovirgula ligni 


81 


0.000248 


0.00132 


Ruegeria atlantica 


8 


0 


0.001062 


Methylosinus sporium 


6 


0 


0.001033 


Methylovirgula ligni 


12 


0 


0.001029 


Pseudorhodoferax soli 


86 


0.000282 


0.00131 


Methylobacterium rhodesianum 


32 


0.000125 


0.001139 


Methylobacterium zatmanii 


22 


0.00001 1 


0.000871 


Flavobocterium psychrophilum 


107 


0.000082 


0.000932 


Rhodobacter capsulotus 


43 


0.000258 


0.001098 


Flavobocterium aquotile 


94 


0.000073 


0.000886 


Rhodobacter sphoeroides 


189 


0.000043 


0.000809 


Methylomicrobium album 


78 


0.000047 


0.000761 


Rhodobacter capsulotus 


91 


0.000226 


0.000929 


Coulobacter leidyia 


46 


0.00002 


0.000712 


Flavobocterium psychrolimnoe 


101 


0.00042 


0.001111 


Brevundimonas voriobilis 


154 


0.00005 


0.000728 


Methylobacterium odhaesivum 


82 


0.000156 


0.000815 


Amaricoccus kaplicensis 


54 


0.000176 


0.000782 


Methylopilo copsulata 


19 


0.000154 


0.000756 


Brevundimonas voriobilis 


71 


0 


0.00054 


Sphingomonas faeni 


182 


0.000046 


0.000579 


Flavobocterium aquotile 


90 


0.000282 


0.000812 


Methylobacterium odhaesivum 


44 


0.000094 


0.000609 


Flavobocterium psychrophilum 


143 


0 


0.000513 


Dokdonella koreensis 


112 


0 


0.00051 1 


Sphingomonas sanguinis 
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Table 3. Cont. 





N° OTU 


mean abundance in group 1 


mean abundance in group 2 


Best hit on Blast 


97 


0 


0.00049 


Hyphomicrobium facile subsp. Tolerans 


105 


0 


0.000359 


Paracoccus haeundaensis 


137 


0 


0.000339 


Acinetobacter johnsonii 


75 


0.000558 


0 


Pseudochrobactrum kiredjianiae 


50 


0.889343 


0.677999 


Methylobacterium rhodesianum 



Differentially abundancy were calculated by using Metastats (using a FDR correction, p<0.01). 
doi:1 0.1 371 /journal.pone.01 02649.t003 



between Melhylobaclerium and Flavobacterium was further 
confirmed by the co-occurrence analysis, which identified poten- 
tial taxonomic interactions between genera. Two significant 
correlations between genera were detected: a positive correlation 
between Maritimibacter and Micrococcus (rj = 0.69, p-value < 
0.01) and a negative correlation between Methylobacterium and 
Flavobacterium (ry=— 0.63, p-value <0.01). Furthermore, vari- 
ance analysis of OTU 36' abundance showed that it was 
influenced by the Fulton index (Shapiro test for Normality: p- 
value =0.07, Linear model: p-value = 0.02741, F= 5.0388, 
Rsq= —0.1496). A negative relationship between Fulton index 
and the OTU 36 was observed (Figure 6). 

In the F2 fish progeny, three significant QTLs (at the genome- 
wide level) were identified on two linkage groups (LG 1 1 and LG 
16). One major QTL per strain was detected respectively for 
Lysobacter, Rheinheimera and Methylobacterium counts (LOD 
score = 9.89, 7.46 and 3.48 respectively). For each of these traits, 
the total PVE (percent variance explained) of the QTL were 
estimated to 17.01%, 31.05% and 41.1% for Methylobacterium, 
Rheinheimera and Lysobacter respectively. The most probable 
positions of these QTL, their respective 95% CIs, the closest linked 
molecular markers (one per QTL) as well as additive and 
dominance effects are presented in Table 6. 

Discussion 

Inter-individual variations in host microbiota has been well 
documented (e.g. [1]. Such variations occur even in hosts with 
identical genetic background, as observed in monozygotic twins 
[3] indicating that both genetic and environmental conditions play 
a role in the modulation of host microbiota. The gold standard 
forward genetics technique to identify areas of the genome that 
relate to certain phenotypes is to make a cross between genetically 
divergent individuals [32]. In this study, we focused on the 
structural variation of the skin microbiota (i.e abundance of each 
bacterial genus) of an F2 intercross originally generated from two 
genetically contrasted strains of brook charr. To our knowledge, 
this is the first study that identified host genomic regions associated 
with the abundance of specific microbiota strains in a non-model 
vertebrate. 

The deep taxonomic analysis of the fish skin microbiota 
indicates that two phyla dominate: Proleobacteria (88.7%), 
followed by Bacteroidetes (9.7%). Those two phyla are mainly 
represented by a single OTU each; 50 and 36 respectively. OTU 
50 is classified as M. rhodesianum. OTU 36 is classified as F. 



psychrophilum, the causative agent of the cold water disease, a 
pervasive infectious disease in farmed fish [33,34]. This disease 
especially affects salmonids at early life-stages, and it is well 
documented that both stressful conditions, and injuries facilitate 
infection triggering [33,35]. According to the dendrogram 
(Figure 2), each individual's bacterial communities clustered into 
three groups. Similar clustering was also found with the second 
normalization method (subsampling of the same number of 
sequences for each samples). The same result was obtained with 
PCoA analysis and furthermore the Pearson correlation indicates 
that the genus Flavobacterium and Melhylobaclerium are nega- 
tively correlated. In the PCoA, two species (F. psychrophilum and 
M. rhodesianum) are correlated in opposite ways with the axis 1 
which discriminate the three groups. We assume that two or more 
species which seem to be mutually exclusive are involved in a 
negative relationship [36]. All together those results indicate that 
the antagonistic relationship between those two dominant species 
is very likely playing a key role in shaping the structure of the 
individual's microbiota. 

Because all fish progenies were reared in the same tank, and 
under the same environmental conditions since birth, environment 
had likely little influence on the phenotypic variation observed in 
this study. Yet, the individual's microbiota either clustered in one 
of two closely related groups, or formed divergent outliers. Sex- 
specific variations had little influence on individual clustering since 
both males and females were present in each defined group and 
there was no correlation between sex and groups. The difference 
between male and female is based on 22 OTUs which are less 
abundant than 0. 1 % of the community and are then classified as 
part of the rare biosphere. Previous studies observed similar 
individual variations of the microbiota in human or mouse 
[1,2,32]. 

Three OTUs classified into two species were responsible for the 
inter-individual differentiation; M. rhodesianum and F. psychro- 
philum. Interestingly, Flavobacterium and Methylobacterium, 
which are the most abundant genera, negatively co-occurred in 
the samples, indicating that the relationship between those two 
genera is based on competition. Furthermore, a strong negative 
correlation was found between the OTU 36 and 50, thus 
supporting a strong antagonistic relationship between the two 
species, F. psychrophilum and M. rhodesianum. 

M. rhodesianum produces poly-P-hydroxybutyrate, a polymer 
of short-chain fatty acid, known to inhibit the growth of pathogens 
like enterobacteria and Vibrio sp. [37,38,39,40,41,42]. Taken 
together, this suggests that a reduction of M. rhodesianum 
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Table 4. Detailed information related to the OTU differentially abundant between individuals belonging to group 1 and group 3. 



OTU 


mean abundance in group 1 


mean abundance in group 3 


Best hit on Blast 


36 


0.00484 


0.15627 


Flavobacterium psychrophilum 


52 


0.000331 


0.052599 


Sphingopyxis olaskensis 


51 


0.002582 


0.045797 


Sphingomonas paucimobilis 


119 


0.001237 


0.034255 


Pseudomonas peli 


125 


0.001567 


0.030681 


Acinetobacter Iwoffii 


120 


0.002302 


0.030821 


Rheinheimera texasensis 


63 


0.00071 


0.027857 


Paracoccus haeundaensis 


122 


0.000778 


0.018741 


Pseudomonas cedrina subsp. Fulgido 


61 


0.003372 


0.02093 


Porphyrobocter dokdonensis 


21 


0.0001 17 


0.01 1826 


Flavobacterium aquatile 


127 


0.000218 


0.00999 


Acinetobacter haemolyticus 


55 


0 


0.008105 


Bradyrhizobium pachyrhizi 


54 


0.000176 


0.007992 


Methylopila capsulata 


162 


0.0021 12 


0.00961 1 


Knoellia sinensis 


126 


0.000969 


0.008421 


Pseudomonas xanthomarina 


67 


0.000518 


0.007951 


Rhodobacter btasticus 


69 


0.000531 


0.007849 


Pseudorhodobacter ferrugineus 


81 


0.000248 


0.006693 


Ruegeria atlantica 


134 


0.000245 


0.00635 


Acinetobacter Iwoffii 


143 


0 


0.005914 


Dokdonella koreensis 


129 


0 


0.005644 


Pseudomonas hibiscicola 


124 


0 


0.004908 


Nevskia soli 


40 


0.000014 


0.004576 


Flavobacterium aquatile 


38 


0.000063 


0.004547 


Sejongia jeonii 


68 


0.0014 


0.005784 


Methylocella palustris 


71 


0 


0.004232 


Sphingomonas faeni 


79 


0.000259 


0.004408 


Rhodobacter capsulatus 


78 


0.000047 


0.003657 


Rhodobacter capsulatus 


37 


0.000471 


0.004035 


Flavobacterium terrigena 


141 


0.000156 


0.003673 


Acinetobacter johnsonii 


156 


0.000033 


0.003487 


Sphingopyxis alaskensis 


39 


0 


0.00318 


Flavobacterium psychrophilum 


24 


0.000095 


0.00295 


Pseudomonas peli 


106 


0.000106 


0.002945 


Bosea vestrisii 


137 


0 


0.002615 


Acinetobacter johnsonii 


101 


0.00042 


0.002997 


Brevundimonas variabilis 


97 


0 


0.002498 


Hyphomicrobium facile subsp. Tolerans 


139 


0.000153 


0.002639 


Lysobacter capsici 


157 


0.00006 


0.002488 


Sphingomonas sanguinis 


116 


0 


0.002339 


Rhodobacter btasticus 


177 


0.000289 


0.00255 


Aquabacterium commune 


25 


0.000069 


0.002236 


Acinetobacter Iwoffii 


182 


0.000046 


0.002172 


Flavobacterium aquatile 


47 


0 


0.002082 


Flavobacterium psychrophilum 


84 


0.000079 


0.00214 


Pseudorhodobacter ferrugineus 


164 


0.000219 


0.002164 


Rhodococcus fascians 


105 


0 


0.001788 


Paracoccus haeundaensis 


43 


0.000258 


0.002022 


Flavobacterium aquatile 


180 


0.000017 


0.001709 


Acidovorax defluvii 


92 


0.000147 


0.001809 


Sphingopyxis witflariensis 
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Table 4. Cont. 



OTU 


mean abundance in group 1 


mean abundance in group 3 


Best hit on Blast 


4 


0 


0.00161 


Methylosinus sporium 


145 


0 


0.001417 


Acinetobacter johnsonii 


5 


0 


0.001415 


Mognetospirillum mognetotacticum 


88 


0.000036 


0.00143 


Amaricoccus veronensis 


3 


0 


0.001358 


Rhodobacter changlensis 


138 


0.000014 


0.001365 


Acinetobacter junii 


70 


0.000831 


0.002052 


Methylobacterium salsuginis 


184 


0 


0.001 156 


Flavobocterium aquotile 


22 


0.00001 1 


0.001 144 


Flavobocterium psychrophilum 


112 


0 


0.001 106 


Sphingomonas sanguinis 


118 


0.000062 


0.001 143 


Haematobacter missouriensis 


93 


0.000029 


0.001093 


Roseomonas mucosa 


152 


0.000028 


0.001009 


Rhodobacter blasticus 


35 


0.00005 


0.000936 


Sphingomonas sanguinis 


148 


0.000033 


0.000915 


Psychrobacter arcticus 


172 


0 


0.000864 


Sphingomonas sanguinis 


146 


0 


0.000832 


Pseudomonas mohnii 


94 


0.000073 


0.000865 


Rhodobacter sphaeroides 


140 


0 


0.000775 


Psychrobacter alimentarius 


160 


0.000061 


0.000823 


Sphingomonas faeni 


83 


0.000128 


0.000877 


Paracoccus aminovorans 


42 


0.000323 


0.001071 


Flavobocterium psychrophilum 


186 


0 


0.00071 1 


Flavobocterium aquotile 


44 


0.000094 


0.000796 


Flavobocterium psychrophilum 


117 


0.000029 


0.000731 


Sphingopyxis taejonensis 


107 


0.000082 


0.000766 


Rhodobacter capsuiotus 


72 


0.000432 


0.001098 


Beijerinckia mobilis 


104 


0.000028 


0.000693 


Sphingopyxis witflariensis 


113 


0 


0.000663 


Stella humosa 


115 


0 


0.000622 


Loktanella salsilacus 


108 


0.000051 


0.000647 


Novosphingobium panipotense 


26 


0 


0.000559 


Acinetobacter johnsonii 


48 


0 


0.00052 


Flavobocterium psychrophilum 


166 


0 


0.000508 


Rhodococcus fascians 


65 


0.00064 


0.001 148 


Methylobacterium organophilum 


103 


0.000215 


0.000655 


Sphingomonas ursincolo 


144 


0.00007 


0.000462 


Acinetobacter Iwoffii 


191 


0.00006 


0.000406 


Beijerinckia derxii subsp. Venezuelae 


19 


0.000154 


0.000498 


Brevundimonas variabilis 


178 


0.000043 


0.000378 


Pseudorhodoferax soli 


23 


0.000017 


0.000348 


Flavobocterium psychrophilum 


192 


0 


0.000328 


Methylovirgula ligni 


49 


0.000058 


0.000354 


Flavobocterium aquotile 


8 


0 


0.00029 


Methylosinus sporium 


9 


0.000062 


0.000245 


Sphingomonas wittichii 


169 


0.000207 


0.000287 


Corynebocterium tuberculostearicum 


111 


0.000416 


0.000202 


Mycoplana bullata 


7 


0.000343 


0 


Singularimonas variicoloris 


175 


0.000348 


0 


Simplicispira limi 
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Table 5. Detailed information related to the OTU differentially abundant between individuals belonging to group 2 and group 3. 



N° OTU 


mean abundance in group 2 


mean abundance in group 3 


Dcsi nil on Blast 


63 


0.002374 


0.027857 


Paracoccus hoeundaensis 


55 


0 


0.008105 


Bradyrhizobium pochyrhizi 


54 


0.000782 


0.007992 


Methylopila capsulata 


134 


0.000155 


0.00635 


Acinetobacter Iwoffii 


129 


0 


0.005644 


Pseudomonas hibiscicola 


143 


0.000513 


0.005914 


Dokdonella koreensis 


124 


0.000109 


0.004908 


Nevskia soli 


40 


0.000306 


0.004576 


Flovobacterium oquatile 


123 


0.000266 


0.004124 


Acinetobacter haemolyticus 


71 


0.00054 


0.004232 


Sphingomonas faeni 


38 


0.001334 


0.004547 


Sejongia jeonii 


156 


0.000306 


0.003487 


Sphingopyxis alaskensis 


39 


0.000144 


0.00318 


Flavobacterium psychrophilum 


78 


0.000761 


0.003657 


Rhodobacter capsulatus 


137 


0.000339 


0.002615 


Acinetobacter johnsonii 


47 


0.000087 


0.002082 


Flavobacterium psychrophilum 


164 


0.000236 


0.002164 


Rhodococcus fascians 


4 


0 


0.00161 


Methylosinus sporium 


182 


0.0005/9 


0.0021 12. 


Flavobacterium aquatile 


3 


0 


0.001358 


Rhodobacter changlensis 


184 


0 


0.001156 


Flavobacterium aquatile 


42 


0.000252 


0.001071 


Flavobacterium psychrophilum 


12 


0.001029 


0 


Pseudorhodoferax soli 


6 


0.001033 


0 


Methylovirgula ligni 


27 


0.001814 


0.000082 


Sphingopyxis chilensis 


50 


0.677999 


0.227483 


Methylobacterium rhodesianum 



Differentially abundancy were calculated by using Metastats (using a FDR correction, p<0.01). 
doi:1 0.1 371 /journal.pone.01 02649.t005 



abundance allows colonization or over-growth by other bacteria, 
some of which are pathogenic. 

Evidently, a change in skin microbiota taxonomic structure 
favors opportunist bacteria and especially opportunistic pathogens 
like F. psychrophilum, Acinetobacter haemolyticus, Acinetobacter 
johnsonii and Acinetobacter junii, which have a higher abundance 
in group 2 compared with group 1 Such a disturbance in the 
microbiome homeostasis is called dysbiosis, and its occurrence 
enlightens the importance of the function of M. rhodesianum in 
controlling the balance between both other commensal bacteria 
and opportunistic pathogens. Furthermore, Flavobacterium psy- 
chrophilum was negatively correlated to both M. rhodesianum 
abundance and Fulton index. The Fulton index is a condition 
factor used as a proxy for fish health status. Therefore, it suggests 
that skin microbiota from the weakest host individuals were unable 
to prevent colonization by the opportunistic pathogen F. 
psychrophilum. Furthermore, skin microbiota taxonomic structure 
varies significantly across the F2 progeny. This result added to the 
fact that i the genetics of the host is the only variable in the 
experiment, and ii the finding of three QTLs associated to the 
abundance of bacterial genera with high PVE, suggest that host 



genotype influences the abundance of commensal strains e.g. 
Methylobacterium, which will regulate the abundance of F. 
psychrophilum. 

We found three QTL associated with the abundance of three 
genera: Lysobacter, Rheinheimera and Methylobacterium, all of 
which were observed to provide antimicrobial compounds 
[39,40,41,42,43,44,45,46,47]. These finding strongly suggests that 
host genotype influences abundance of specific commensal strains, 
and possibly targets their recruitment. The most compelling 
evidence concerns the QTL associated with Methylobacterium 
(PVE= 17.01%): as previously described above, Methylobacterium 
is influential on the structure and the homeostasis of the 
microbiota. Its abundance is inversely correlated with those of 
the pathogen F. psychrophilum. The targeted recruitment of M. 
rhodesianum mediated by the host genotype is therefore a strategy 
to prevent pathogen growth via harnessing antagonistic relation- 
ship towards resources use [48]. To our knowledge, this is only the 
second study that identified QTL associated with microbial 
variation among individuals. A study on murine gut microbiota 
showed similar results, in which McKnite et al. (2012) found 5 
QTL located on four chromosomes influencing the variation of 
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Figure 5. PCoA analysis of the microbiome for all 86 F 2 individuals based on the Yue & Clayton measure of dissimilarity (Thetayc). 

Black circles represent individuals belonging to the group 1, red circles represent individuals belonging to the group 2 and green circles represent 
individuals belonging to the group 3. Arrows represent genus, which are significantly correlated with the axis. The first and second axes represented 
24.5% and 7.7% of the variation respectively. The R-squared between the original distance matrix and the distance between the points in 2D PCoA 
space was 0.87. 

doi:1 0.1 371 /journal.pone.01 02649.g005 



different taxa. The variance explained by their QTL is in the same 
range of our (20% to 27% of the variance explained for McKnite 
el al. (2012)). We also found QTLs with a major effect on the 
variance of genus abundances (PVE ranging for 17.02 to 41.1%). 
Linkage analysis on genus abundance data strongly evidenced the 
influence of host genetics on the modulation and/ or recruitment 
for those genera. Therefore, brook charr immunity involves both 
specific commensal strains recruitment and nuclear gene expres- 
sion, as previously evidenced in [49], those are under the control 
of the host genotype. 

The study of the genetic architecture underlying the regulation 
of bacterial abundance further highlights the coevolutionary 



pattern between host, commensals, and pathogens. However, 
identifying genes located in the genomic regions (QTL) linked to 
the abundance of microbial partners will be far more challenging 
as the genome of brook charr is currendy not fully sequenced and 
poorly annotated compared to the mouse genome. The markers 
surrounding the QTL regions may define regions to be deeply 
sequenced in future work to identify potential underlying genes 
and their associated functions. Evidently, those markers will be 
invaluable to conduct highly innovative genetic breeding programs 
targeting disease resistant host strains via the recruitment of highly 
resilient microbiota. 
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Table 6. Descriptive statistics, including the LOD score, the position, 95% CI, PVE (%), the associated P-value of each QTL linked to 
every phenotype related to bacteria counts trait (LOD, Log qo of the odd ratio; 95% CI, 95% confidence interval; PVE, percent 
variance explained) [23]. 





Phenotype 


Linkage Group 


Pos (cM) 


95% CI (cM) 


LOD 


p value (F) 


PVE (in %) 


Nearest marker 


Lysobacter 


11 


62 


56.2-68.2 


9.89 


0.000 


41.1 


sf003455 


Rheinheimera 


16 


42.8 


39.3-45.5 


7.46 


0.000 


31.05 


SFO-D91 


Methylobacterium 


16 


39 


34.5-43.5 


3.48 


0.001 


17.02 


SFO-D91 
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Supporting Information 

Figure SI Rarefaction curves for each group defined by 
the dendrogram analysis based upon ThetaYC index 
(see figure 2). Each group curve reaches a plateau thus 
indicating the depth of sequencing is sufficient. 
(EPS) 

Figure S2 Alpha diversity of each group normalized by 
subsampling of the same number of reads per sample. 

Each group was defined by the dendrogram analysis based upon 
ThetaYC index performed on the normalized dataset. Alpha- 
diversity was calculated by the non-parametric index of Shannon. 
Statistical differences (represented by the asterisk) were calculated 
with a wilcoxon test with a correction for multiple testing (Holm, 
p-values <1.10 4 ). 
(EPS) 

Figure S3 PCoA analysis of the microbiome for all 86 F 2 
individuals. This PCoA was performed on a normalized dataset 
with the second method of normalization (subsampling of the same 
number of reads per sample). Black circles represent individuals 
belonging to the group 1, red circles represent individuals 
belonging to the group 2 and green circles represent individuals 
belonging to the group 3. Arrows represent genera, which are 

References 

1. Costcllo EK, Laubcr CL, Hamady M, Fierer N, Gordon JI, ct al. (2009) 
Bacterial community variation in human body habitats across space and time. 
Science 326: 1694-1697. 

2. Fierer N, Laubcr CL, Zhou N, McDonald D, Costcllo EK, ct al. (2010) From 
the cover: Forensic identification using skin bacterial communities. Proceedings 
of the National Academy of Sciences 107: 6477—6481. 

3. Turnbaugh PJ, Quince C, Faith JJ, McHardy AC, Yatsunenko T, et al. (2010) 
Organismal, genetic, and transcriptional variation in the deeply sequenced gut 
microbiomes of identical twins. Proceedings of the National Academy of 
Sciences 107: 7503-7508. 

4. Turnbaugh PJ, Ley RE, Hamady M, Fraser-Liggett CM, Knight R, et al. (2007) 
The human microbiome project. Nature 449: 804—810. 

5. Turnbaugh PJ, Henrissat B, Gordon JI (2010) Viewing the human microbiome 
through three-dimensional glasses: integrating structural and functional studies 
to better define the properties of myriad carbohydrate-active enzymes. Acta 
Crystallographica Section F 66: 1261—1264. 

6. Turnbaugh PJ, Gordon JI (2009) The core gut microbiome, energy balance and 
obesity. The Journal of Physiology 587: 4153-4158. 

7. Giongo A, Gano KA, Crabb DB, Mukherjcc N, Novelo LL, ct al. (2011) 
Toward defining the autoimmune microbiome for type 1 diabetes. ISMEJ 5: 
82-91. 

8. Peterson DA, McNulty NP, Guruge JL, Gordon JI (2007) IgA Response to 
Symbiotic Bacteria as a Mediator of Gut Homeostasis. Cell Host & Microbe 2: 
328-339. 

9. Bouskra D, Brezillon C, Berard M, Wcrts C, Varona R, ct al. (2008) Lymphoid 
tissue genesis induced by commensals through NODI regulates intestinal 
homeostasis. Nature 456: 507-510. 

10. Macphcrson AJ, Geuking MB, McCoy KD (201 1) Immunoglobulin A: a bridge 
between innate and adaptive immunity. Current Opinion in Gastroenterology 
27: 529-533 510.1097/MOG.1090bl013e32834bb32805 

1 1 . Boutin S, Bernatchcz L, Audet C, Derome N (2012) Antagonistic effect of 
indigenous skin bacteria of brook charr (Salvelinus fontinalis) against 
Flavobacterium columnare and Fpsychrophilum. Veterinary Microbiology 155: 
355-361. 

12. Balcazar J, de Bias I, Ruiz-Zarzuela I, Vendrcll D, Muzquiz J (2004) Probiotics: 
a tool for the future of fish and shellfish health management. Journal of 
Aquaculturc in the Tropics 19: 239-242. 

13. Backhed F (2005) Host-bacterial mutualism in the human intestine. Science 307: 
1915-1920. 

14. Hooper LV, Gordon JI (2001) Commensal host-bacterial relationships in the 
gut. Science 292: 1115-1118. 

15. O'Hara AM, Shanahan F (2006) The gut flora as a forgotten organ. EMBO 
Reports 7: 688-693. 

16. Stecher B, Hardt W-D (2008) The role of microbiota in infectious disease. 
Trends in Microbiology 16: 107-114. 

17. Loiscau C, Richard M, Gamier S, Chastcl O, Julliard R, ct al. (2009) 
Diversifying selection on MHC class I in the house sparrow (Passer domesticus). 
Molecular Ecology 18: 1331-1340. 



significandy correlated with the axis. The first and second axes 
represented 24.5% and 7.8% of the variation respectively. The R- 
squared between the original distance matrix and the distance 
between the points in 2D PCoA space was 0.88. 
(EPS) 

Table SI Differentiation between the three groups of 
OTU calculated with an analysis of similarity (ANOSIM, 
statistic R). The q-values represent p-values corrected with 
Bonferroni and considered as significant when p<0.05. Note: The 
analysis was performed on a dataset normalized by two methods: 
zscore and subsampling. Both analyses gave the same results, each 
group being significandy different from the two others. 
(DOC) 

Acknowledgments 

We thank Martin Llewellyn, Scott A. Pavey, the associate editor of Plos 
One, and two anonymous reviewers for critical reading of the manuscript. 

Author Contributions 

Conceived and designed the experiments: ND CA LB SB CS. Performed 
the experiments: SB GS. Analyzed the data: SB CS. Contributed reagents/ 
materials/analysis tools: ND CA LB. Wrote the paper: ND CA LB SB CS. 



18. Krawczyk M, Reith W (2006) Regulation of MHC class II expression, a unique 
regulatory system identified by the study of a primary immunodeficiency disease. 
Tissue Antigens 67: 183—197. 

19. Dionnc M, Miller KM, Dodson JJ, Bernatchcz L (2009) MHC standing genetic 
variation and pathogen resistance in wild Atlantic salmon. Philosophical 
Transactions of the Royal Society B-Biological Sciences 364: 1555-1565. 

20. Alexander C, Rietschel ET (2001) Invited review: Bacterial lipopolysaccharidcs 
and innate immunity. Journal of Endotoxin Research 7: 167—202. 

21. Boutin S, Audet C, Derome N (2013) Probiotic treatment by indigenous bacteria 
decreases mortality without disturbing the natural microbiota of Salvclinus 
fontinalis. Canadian Journal of Microbiology In press. 

22. Boutin S, Bernatchcz L, Audet C, Derome N (2013) Network Analysis 
Highlights Complex Interactions between Pathogen, Host and Commensal 
Microbiota. PLoS ONE 8: c84772. 

23. Sauvage C, Vaguer M, Derome N, Audet C, Bernatchcz L (2012) Coding gene 
single nucleotide polymorphism mapping and quantitative trait loci detection for 
physiological reproductive traits in brook charr, Salvelinus fontinalis. G3: 
Genes | Genomes | Genetics 2: 379—392. 

24. Livia L, Antonella P, Hovirag L, Mauro N, Panara F (2006) A nondestructive, 
rapid, reliable and inexpensive method to sample, store and extract high-quality 
DNA from fish body mucus and buccal cells. Molecular Ecology Notes 6: 257— 
260. 

25. Turner S, Prycr KM, Miao VPW, Palmer JD (1999) Investigating deep 
phylogenetic relationships among C,\anobacteria and plastids by small subunit 
rRNA sequence analysis. Journal of Eukaryotic Microbiology 46: 327—338. 

26. Marchesi JR, Sato T, Wcightman AJ, Martin TA, FryJC, et al. (1998) Design 
and evaluation of useful bacterium-specific PGR primers that amplify genes 
coding for bacterial 16S rRNA. Applied and Environmental Microbiology 64: 
795-799. 

27. Schloss P, Westcott S, Ryabin T, Hall J, Hartmann M, ct al. (2009) Introducing 
mothur: open-source, platform-independent, community-supported software for 
describing and comparing microbial communities. Applied and Environmental 
Microbiology 75: 7537-7541. 

28. Siddharth J, Holway N, Parkinson SJ (2013) A Western Diet Ecological Module 
Identified from the 'Humanized' Mouse Microbiota Predicts Diet in Adults and 
Formula Feeding in Children. PLoS ONE 8: c83689. 

29. Takcshita T, Suzuki N, Nakano Y, Yasui M, Yoncda M, ct al. (2012) 
Discrimination of the oral microbiota associated with high hydrogen sulfide and 
methyl mercaptan production. Sci Rep 2. 

30. Horner-Devine MC, Lage M, Hughes JB, Bohannan BJM (2004) A taxa-area 
relationship for bacteria. Nature 432: 750-753. 

31. White JR, Nagarajan N, Pop M (2009) Statistical methods for detecting 
differentially abundant features in clinical metagenomic samples. PLoS 
Computational Biology 5: clOOO.352. 

32. McKnite AM, Perez-Munoz ME, Lu L, Williams EG, Brewer S, ct al. (2012) 
Murine gut microbiota is defined by host genetics and modulates variation of 
metabolic traits. PLoS ONE 7: c.39191. 

33. Starliper CE (2011) Bacterial eoldwater disease of fishes caused by Flavobac- 
terium psychrophilum. Journal of Advanced Research 2: 97-108. 



PLOS ONE | www.plosone.org 



16 



July 2014 | Volume 9 | Issue 7 | e102649 



Genotypes Control Fish Skin Microbiota 



34. Wiklund T, Madsen L, Bruun MS, Dalsgaard I (2000) Detection of 
Flavobacterium psychrophilum from fish tissue and water samples by PGR 
amplification. Journal of Applied Microbiology 88: 299-307. 

35. Madetoja J, Nyman P, Wiklund T (2000) Flavobacterium psychrophilum, 
invasion into and shedding by rainbow trout Oncorhynchus mykiss. Diseases of 
Aquatic Organisms 43: 27-38. 

36. Faust K, RacsJ (2012) Microbial interactions: from networks to models. Nat Rev 
Micro 10: 538-550. 

37. Defoirdt T, Halet D, Vervaercn H, Boon N, Van dc Wide T, et al. (2007) The 
bacterial storage compound poly-p-hydroxybutyratc protects Artemia francis- 
cana from pathogenic Vibrio campbellii. Environmental Microbiology 9: 445- 
452. 

38. Halct D, Dcfbirdt T, Van Damme P, Vervaercn H, Forrez I, et al. (2007) Poly- 
P-hydroxybutyrate-aceumulating bacteria protect gnotobiotic Artemia francis- 
cana from pathogenic Vibrio campbellii. FEMS Microbiology Ecology 60: 363- 
369. 

39. Yellore Desai (1998) Production of poly-3 -hydro xybutyrate from lactose and 
whey by Methylobactermm sp. ZP24. Letters in Applied Microbiology 26: 391— 
394. 

40. Mothes G, Ackermann J-U, Babel W (1998) Regulation of poly(p-hydroxybu- 
tyrate) synthesis in Methylobacterium rkodesianum MB 126 growing on methanol 
or fructose. Archives of Microbiology 169: 360-363. 

41. Bormann EJ, Roth M (1999) The production of polyhydroxybutyrate by 
Methylobacterium rhodesianum and Ralstonia eutropha in media containing 
glycerol and casein hydrolysates. Biotechnology Letters 21: 1059-1063. 



42. Bormann EJ, Roth M, LinR W, Lcipncr M (1999) Selection of capsule deficient 
strains of Methyl 'ohactenum rhodesianum producing polyhydroxybutyrate. 
Biotechnology Techniques 13: 539-544. 

43. Balachandran C, Duraipandiyan V, Ignaeimuthu S (2012) Gytotoxie (A549) and 
antimicrobial effects of Methylobacterium sp. isolate (ERI-135) from Nilgiris 
forest soil, India. Asian Pacific Journal of Tropical Biomedicinc 2: 712—716. 

44. Korotkova N, Chistoserdova L, Lidstrom ME (2002) Poly-p-Hydroxybutyratc 
Biosynthesis in the Facultative Mcthylotroph Methylobacterium cxtorquens 
AMI: Identification and Mutation of gapll, gap20, and phaR. Journal of 
Bacteriology 184: 6174-6181. 

45. Chen W-M, Lin C-Y, Sheu S-Y (2010) Investigating antimicrobial activity in 
Rheinheimera sp. due to hydrogen peroxide generated by 1-lysinc oxidase 
activity. Enzyme and Microbial Technology 46: 487^£93. 

46. Chen W-M, Lin C-Y, Young C-C, Sheu S-Y (2010) Rheinheimera aquatica sp. 
nov., antimicrobial activity- producing bacterium isolated from freshwater 
culture pond. Seoul, COREE, REPUBLIQUE DE: Korean Society for Applied 
Microbiology. 7 p. 

47. Reichenbach H (2006) The genus Lysobacter. In: Dworkin M, Falkow S, 
Rosenberg E, Schlcifcr K-H, Stackebrandt E, editors. The Prokaryotcs: Springer 
New York. pp. 939-957. 

48. Dillon R, Charnley K (2002) Mutualism between the desert locust Schistocerca 
gregaria and its gut microbiota. Research in Microbiology 153: 503-509. 

49. Langevin C, Blanco M, Martin SAM, Jouncau L, Bernardet J-F, et al. (2012) 
Transcriptional responses of resistant and susceptible fish clones to the bacterial 
pathogen Flavobacterium psychrophilum. PLoS ONE 7: c39126. 



PLOS ONE | www.plosone.org 



17 



July 2014 | Volume 9 | Issue 7 | e102649 



